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Abstract. We develop a method to study the anisotropy of a cosmic-ray angular distribution, using both the right ascension and 
the declination of the arrival directions. It generalises the full-sky coverage method of Sommers (2001) to partial-sky coverage 
experiments. When the angular distribution consists of a dipolar modulation of an otherwise isotropic flux, the method allows 
one to reconstruct the dipole amplitude and the dipole orientation in 3D space. We analyse in detail the statistical properties 
of the method, introducing the concept of reconstruction power, and show that it is generally more powerful than the standard 
Rayleigh analysis in right ascension. We clarify the link between the traditionally-used first harmonic amplitude and the true, 
physical dipole amplitude, and we investigate the variation of the reconstruction powers as a function of the dipole orientation. 
We illustrate the method by computing the amplitude and angular reconstruction powers of the Pierre Auger Observatory, 
with the Southern site alone and with both Southern and Northern sites. In this particular case, we find that the building of a 
similar site in the Northern hemisphere would decrease the time needed for the method to reveal a significant departure from 
an isotropic cosmic-ray distribution by a factor of about eight. 
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1. Introduction 

Although cosmic-rays (CRs) have been known for almost one 
century, their origin remains uncertain, mostly because their 
trajectories are bent by Galactic magnetic fields and they do 
not individually point back to their sources. Moreover, since 
these fields are chaotic on scales ranging at least from 10 8 cm 
to 10 20 cm (Armstrong et al. 1995), the transport of CRs is dif- 
fusive up to high energies, which tends to make their angular 
distribution isotropic. Therefore, even collectively, the CR ar- 
rival directions hold virtually no information about the source 
distribution in space. 

However, as the energy of the CRs increases, anisotropies 
can appear either because the diffusive approximation does not 
hold anymore, or because the diffusion coefficient becomes 
large enough to reveal intrinsic inhomogeneities in the source 
distribution. Specifically, even if the diffusive regime holds, the 
density of CR sources in the Galaxy is believed to be larger in 
the inner regions than in the outer ones, and this can cause a 
slightly higher CR flux coming from the Galactic centre (GC) 
than from the anti-center. Likewise, the global CR streaming 
away from the Galactic plane (towards the halo) can be a source 
of measurable anisotropy. However, the detailed angular dis- 
tribution of CRs is quite hard to predict, even if we assume 
a definite source distribution, because it also depends on the 
propagation conditions, which are related to both large scale 
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and small scale magnetic field configurations, and on the posi- 
tion of the Earth relative to major magnetic structures, such as 
the local Galactic arm. The presence of numerous nearby su- 
perbubbles, which can break out the Galactic disk and produce 
chimneys (e.g. de Avillez & Berry 2001) along which CRs dif- 
fuse more easily into the halo, can also be responsible for spe- 
cific CR anisotropies. 

Under the simplest assumptions (isotropic diffusion, homo- 
geneous distribution of sources in the Galactic disk), the off- 
centered position of the Earth in the Galaxy (radially and verti- 
cally) leads to a dominant dipole anisotropy with an amplitude 
a proportional to the energy-dependent diffusion coefficient, 
D(E) (e.g. Jones et al. 2001). At higher energy, it is also ex- 
pected that the influence of local structures will become less 
important. A dipole anisotropy, although with a small ampli- 
tude, can also arise because of the relative motion of the solar 
system with respect to the interstellar plasma carrying the dif- 
fusing centers. At ultra-high-energy (UHE), a dipole can also 
appear if CRs propagate in straight lines from exotic sources 
distributed all over the Galactic halo, e.g. related to the dark 
matter (e.g. Berezinsky 2000), or if a roughly diffusive propa- 
gation settles between a dominant source (e.g. the Virgo clus- 
ter) and our Galaxy. Higher order multipoles can also arise nat- 
urally in some models. For instance, a dominant quadrupole in 
the UHECR angular distribution could result from an equato- 
rial excess of sources in supergalactic coordinates (Sommers 
2001). 
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From a general point of view, the characterisation of the 
CR anisotropy provides useful information to constrain the 
GCR diffusion models, notably the effective diffusion coeffi- 
cients, related to the magnetic field structure. Indeed, the level 
of CR anisotropy depends on the diffusion coefficient: in a sim- 
ple model where CR sources are homogeneously distributed 
in a disk of thickness 2h and the CRs are confined in a halo 
of height H, the anisotropy at a distance z above the Galactic 
plane (z < h) is estimated as 5 ^ 3D/cH x z/h (Ptuskin 1997). 
Anisotropy measurements at various energies can thus provide 
crucial information about the energy dependence of the dif- 
fusion coefficient. This information is particularly important 
to constrain the GCR source spectrum, since it sets the re- 
lation between the source power-law index and the observed 
one, through the energy dependent confinement of CRs in the 
Galaxy. 

In Sect. we present a method to derive the parameters of 
an assumed dipole anisotropy in the CR angular distribution, 
following Sommers (2001) in using both the right ascension, 
RA, and the declination, 6, of a discrete set of CRs. This is in 
contrast with the standard "Rayleigh analysis" (i.e. harmonic 
analysis on the circle) used in cosmic-ray physics, which makes 
use of the CR distribution in right ascension only (Linsley 
1975). The method exposed below consists essentially of a gen- 
eralization of Sommers' method (i.e. harmonic analysis on the 
sphere, Sommers 2001) to limited sky coverage. In Sect. [5] we 
study quantitatively the statistical power of the method for the 
reconstruction of the dipole amplitude and orientation, in the 
3-dimensional space (two angular dimensions). We then apply 
this study to the case of the Pierre Auger Observatory (PAO, 
e.g. Auger Collaboration 2004), which will soon obtain the 
largest data set in the crucial energy range between 10 18 and 
5 10 18 eV, i.e. at the end of the GCR component, before the 
transition to extragalactic CRs. We then compare the accuracy 
of the Rayleigh analysis with that of our 3D method. We also 
discuss the increase in the reconstruction power which will be 
provided by the second site of the PAO, to be installed in the 
Northern hemisphere. 

2. Dipole reconstruction procedure 

2.1. Anisotropy in right ascension: the Rayleigh 
analysis 

The standard way to estimate the anisotropy of the CR distribu- 
tion at a given energy consists in performing a Fourier analysis 
of the CR arrival directions and computing the first-harmonic 
amplitude in right ascension, r^. This is obtained from the 
sums of the sine and cosine of the right ascension of each of 
the N events, which is also the azimuthal angle, <p, in equato- 
rial spherical coordinates: 

2 ft 2 N 

a=—^cos<pi and ^sin^;, (1) 

/=i /=i 

as 

r lh = Va 2 + b 2 . (2) 



The direction of the anisotropy is also obtained from the phase 
of the first harmonic in right ascension: 

b 

iffih = arctan -. (3) 
a 

The reconstructed values of the first harmonic amplitude 
and phase converge towards the actual ones when the number 
of CRs tends to infinity. The distribution of reconstructed val- 
ues for limited data sets has been analysed formally by Linsley 
(1975). For sufficiently large numbers of events, N » 4/rj h , 
it is well described by a Gaussian distribution with widths 

cr r = yJTjN and cr^ = ^2/Nr 2 h . For the purpose of compari- 
son with the 2D approach described below, we can reformulate 
this result by defining the significance of the measurement of a 
given anisotropy amplitude, r^, as the "number of sigmas": 

n (r = — = ^ F r lb VA7 (4) 
cr r a/2 

With the above method, the direction of the flux excess is 
only known by its RA coordinate, since no information about 
the declination is available. This is an obvious consequence of 
the initial choice of marginalising with respect to this variable, 
which is usually made for practical reasons. The acceptance 
of a CR experiment usually depends much more on declina- 
tion than on right ascension, because all the points of the sky 
at the same declination are observed in roughly the same con- 
ditions as the Earth rotates around its poles. In the case of ex- 
periments working around 1-100 TeV, the CR anisotropy is so 
small (see Sect. Q that even a limited uncertainty of the order 
of 1 % in the declination dependence of the acceptance would 
completely blur the results. 

At higher energies, however, the anisotropy is expected 
to be higher, and the determination of its direction, in right 
ascension and declination, appears both possible and desir- 
able. In the next subsection, we recall the method proposed by 
Sommers (2001) in the case of a full sky coverage experiment, 
and then we give a generalization for partial-sky coverage. 

2.2. Complete dipole reconstruction with full sky 
coverage 

Let us start by writing the assumed shape of the CR angular 
distribution: 

O(u) = — (1 + »D u), (5) 
4n 

where the differential flux, in direction u, consists of an 
isotropic part, €>o/47r, modulated by a dipolar component, in 
cos(u,D). Here, D is the unit vector pointing in the direction 
of the dipole, and a is the dipole amplitude, relative to the 
monopole: < a < 1. This can be seen as the first order devel- 
opment in spherical harmonics of the CR angular distribution. 
In the following, we shall always assume that the higher order 
terms (quadrupole, etc.) are negligible, so that the flux which 
we want to reconstruct is of the form of Eq. 10. Specifically, 
we want to derive a and D from the data, i.e. three parameters: 
one amplitude, and two angles. 

The idea is to compute some integral quantities and to iden- 
tify their theoretical values with the discrete versions provided 
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by discrete sums over experimental data. In order to reconstruct 
the three parameters of a dipole, we need three quantities and 
an additional one corresponding to the global flux normaliza- 
tion. This is provided by the following moments of order zero 
and one: 



k = J <E>(u)dQ and I = j u<D(u)d£l 



(6) 



It is straightforward to obtain these integrals over the whole 
sky when <I>(u) is given by Eq. 



1 

/ = ® and I = -O x arD. 



(7) 



The discrete version of these integrals, So and S, are ob- 
tained by dividing the sky into a series of pixels {(/, j)}, with 
solid angle Scoij, and changing the continuous integral into a 
sum over all the pixels: 



/(u)dQ 



(U) 



(8) 



Let Nij be the number of CRs observed in the direction of 
the pixel (i, j). Introducing the exposure of a given direction 
on the sky, £(u) (in m 2 s), relevant to a particular CR experi- 
ment, one can write the differential number of events observed 
in direction u, within dQ, as: 



dN . 
dQ 



(u)dO = <t>(u)«5(u)dQ — > N u = Q> u 8 u S<o, 



(9) 



Transforming the integrals in Eq. l|6} according to (|SJl, and 
replacing Su>ij from Eq. one obtains: 



s -Y NiJ 



and I 



(ij) ,J 



(10) 



where Uy is the unit vector in the direction of pixel (i, j). 

The last step consists in changing the sum over all direc- 
tions (or pixels) into a sum over all events. This is done by not- 
ing that when one sums over all the events in the same pixel, 
(i, j), one actually adds up the corresponding value N,j times. 
Therefore, for any quantity F, if F,-j denotes its value in direc- 
tion (;, j), and Fk is its value in the direction of event k, one can 
write: 

(id 

(ij) k 

where the first sum is over all the pixels, and the second is over 
all the events. 

Applying this to Eqs. I ll PL one finally gets the discrete ver- 
sions of Eq. (|6}: 



and 



(12) 



where the sums are over all the events and is the exposure of 
the sky in the direction of event k, namely Ujt, as observed with 
the experiment under consideration. 

The above four discrete sums, So and S, can be straightfor- 
wardly computed from the data (provided that the sky exposure 
is known, which simply derives from the detector's aperture). 



The derivation of the dipole parameters then follow directly 
from the identification with 7o and I, as given by Eqs. 0. The 
dipole amplitude and directions are estimated (given the finite 
set of events available) as: 



3^ and D : 

So 



_S_ 



(13) 



Obviously, the above estimates are all the more accurate 
that the total number of events, N, is large (see Sect. [3}- 



2.3. Generalization to the case of partial sky coverage 

As noted by Sommers (2001), the discrete sums of Eq. fllOi 
are meaningless in regions where the exposure is null. This 
prevents one from using the method in its current form with 
partial-sky coverage data sets. However, the above derivation 
makes it clear that a generalisation to such a case is possible, 
although slightly more complicated from the algebraic point of 
view. The simplicity of the method was due to the fact that the 
matrix giving 7o and I as a function of the dipole parameters, 
a and D, could be straightforwardly inverted to give a and D 
as a function of the integral moments, or their approximations 
So and S, as summed up in Eqs. Jl 3i . Exactly the same pro- 
cedure can be followed in the case of a limited sky coverage 
experiment. 

Since the sums in Eqs. dl 31 are over detected events, it is 
clear that the exposures, in the denominator can never be 
zero, even though there are parts of the sky where the exposure 
is null. To identify So and S with the corresponding moments, 
/() and I, one simply needs to integrate the latter only over the 
part of the sky which is actually observed. Interestingly, the sky 
region covered by any terrestrial observatory (operating for a 
long enough period of time) will always be limited by lines of 
constant declination, i.e. independent of RA. This is an obvious 
consequence of the Earth's rotation around its axis and the Sun. 
In the following, we will thus assume that the part of the sky 
where the exposure is non-zero is contained between declina- 
tion S m i n and <5 max (in equatorial coordinates), corresponding to 
spherical 6 coordinates (as measured from the North Pole) be- 
tween 9 m i n and # max . Note that we do not require the exposure 
to be independent of RA (which is usually not the case, because 
of seasonal effects, e.g. related to temperature variations of the 
acceptance), but only that if a point on the celestial sphere with 
declination 5 can be observed by the detector, then all the points 
at the same declination can also be observed (possibly with a 
different exposure), whatever their right ascension. 

We thus rewrite the zeroth and first order moments of the 
CR intensity as: 



/o 



d9 sin0 \ dy><t>(u) 



JB min JO 
r*8 max r-ln 

and I=| d<9sin<9 dy3U<I>(u), 

Je mi „ Jo 



(14) 
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which can be integrated using Eq. (0 as: 



/o 

h 

h 
I- 



<&od 

-f-(2 + saD 7 ) 
4 

-^-(5 + 2yaD z ) 



(15) 



where we have introduced the difference, sum and product 

d — COS Oroia — COS # max 

s = cos 6> mi n + cos 6> max (16) 
p = cos min x cos 6» max 



and defined 



7 = 



2 

s L - p 



(17) 



Finally, Eqs. dl5> can be inverted to obtain 



aD x 
aD y 
aD 7 = 



y-p 

sl z - 2yl y - 1 

h y-p 

sl z — 2yIo y - 1 
sip ~ 2/ z 
sl z - 2yIo 



(18) 



As in the case of a full sky coverage, the dipole parameters 
can thus be evaluated straightforwardly from the discrete ver- 
sions of /(> and I, which are still given by Eqs. H2\ . The dipole 
direction is reconstructed as the direction of the vector on the 
right hand side of Eqs. (II 81 . while its amplitude is its norm. 

For instance, the reconstructed dipole amplitude, a rec , 
writes: 



ly — 



2s z y 



(2yS - sS z ) 2 



(19) 



which reduces to a 2 ec = 9(S X + S 2 + S z )/Sq, as expected, in the 
case of a full sky coverage experiment (s = 0, d = 2, p = -1, 
7=1/3). 

An obvious advantage of this method with respect to the 
2D Rayleigh method is that it entirely characterises the dipole 
vector, i.e. the direction of the anisotropy both in declination 
and right ascension. One may be worried, however, that this is 
at the expense of reconstruction accuracy, since it is essentially 
the same data set which is used to derive one more (angular) in- 
formation. We show in the next section that this is not the case, 
and that the dipole reconstruction accuracy is in fact generally 
better with the method presented here. 

3. Reconstruction accuracy 

3.1. Monte-Carlo technique 

The power of Eqs. (II 81 for a dipole reconstruction can be eval- 
uated by generating random samples of events, and building the 
sums So and S in Eq. (1121 to reconstruct the underlying dipole 




Fig. 1. Relative exposure of the PAO surface detector (either 
one or two sites) as a function of declination in equatorial co- 
ordinates, assuming efficient detection up to a maximum zenith 
angle of either 60° or 80°. 



blindly (cf. Sommers, 2001). We thus generated artificial data 
sets of various sizes, according to the anisotropic flux given in 
Eq. taking also into account the relative exposure of the 
different parts of the sky as observed by any particular experi- 
ment of our choice. This comes down to drawing randomly N 
directions over the sky, i.e. N pairs of angles (0,, ipi) represent- 
ing the arrival direction of the cosmic -rays, with a probability 
fiyi) = O(u) x £(u), where <£>(u) is given by Eq. (0 and £(u) 
depends on the experiment considered. 

As an illustration, we consider the case of the Pierre Auger 
Observatory (PAO), either with the Southern site only or with 
both the Southern and Northern sites. We first calculate the cor- 
responding exposure function for CRs at sufficiently high en- 
ergy (so that the detection efficiency is saturated, i.e. indepen- 
dent of the CR arrival direction). The effective detection area is 
thus simply the ground area multiplied by cos 8 Z , where 6 Z is the 
local zenith angle. The exposure of a given point in the "equa- 
torial sky" (declination, S, and right ascension, RA) is then ob- 
tained by integrating over the path of that point in the "local 
sky" over the detector, which depends on its latitude on Earth. 
The Southern site is located at A ^ -35.2°, while the Northern 
site is not chosen yet, but should be at a similar (positive) lat- 
itude. The result is shown in Fig.^ f° r different assumptions 
concerning the maximum zenith angle up to which the detec- 
tor is assumed to be working, namely 6* z max = 60° and 80°. 
The PAO is expected to be working in standard reconstruction 
mode up to ~ 70°, which is the value which we adopt for 6> Zimax 
in the calculations below (Auger Collaboration, 2004). 

We then estimate the dipole reconstruction accuracy for 
various values of the dipole amplitude, a, and direction, (6i,(pi), 
and for various sizes of the data set, N. In each case, we simu- 
late a large number of data sets (namely 10 4 ), and compare the 
reconstructed dipole parameters with the input values. 
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Fig. 2. Left: distributions of the relative error on the reconstructed dipole amplitude for various sizes of the data set, as indicated, 
in the case of the PAO detector with Southern site only, and a maximum zenith angle of 70°. Right: bias and dispersion of the 
reconstructed amplitude for a dipole with a = 0.10 pointing towards the Galactic centre (GC), as a function of the number of 
events. 
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Fig. 3. Variation of the dispersion of Sal a, as a function of a, for various sizes of the data set (left) and as a function of N, for 
various values of a. The fit given by formula ( I20I is also drawn. 



3.2. Reconstruction of the dipole amplitude and 
direction 

The distribution of the reconstructed dipole amplitudes are 
shown on Fig.^k for various values of N, in the case of a dipole 
oriented towards the GC with an amplitude of 10% (a - 0.1). 
The plotted quantity is the relative error on the reconstructed 
amplitude, namely 6a /a = (a tme - a KC )/a tme . Its distribution 
is seen to be very close to Gaussian, even for relatively small 
data sets. The bias and dispersion both decrease when the data 
set gets larger. In Fig. |2j), this decrease is shown to follow re- 
spectively a A^ 1 law and a N~ l/2 law at large N, as expected 
from statistics. Note also that the bias is always smaller than 
the dispersion. 



On Fig. |3] we show the dispersion of 6a I a as a function 
of a, for various sizes of the data set, and as a function of N, 
for various values of a. As can be seen, it obeys a general law 
cr(6a/a) oc (a yfRy 1 , which allows us to write the significance 
of a given measurement (number of sigmas, as in Eq.|3: 

n a = [o-(6ala)Y l = K a a^fN, (20) 

where we have introduced the so-called reconstruction power, 
K a , of the detector. This parameter fully characterises the abil- 
ity of the detector under consideration to measure the ampli- 
tude of an underlying dipole anisotropy in the CR flux. In the 
case of a non uniform sky coverage, the reconstruction power 
also depends on the orientation of the dipole, as further dis- 
cussed in the next subsection. For instance, in the case of a 
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dipole oriented in the GC direction, observed by the Southern 
site of the PAO, the reconstruction power is K a = 0.33, as 
shown by the fits on Fig. [2] 

Likewise, we show on Fig.|4t the distribution of the error 
on the reconstructed dipole declination, 5 - n/2-9, for various 
sizes of the data set, and in Fig. 0f), the evolution of the bias 
and dispersion in 5 (similar results are obtained with the other 
angular dimension, namely the right ascension of the dipole 
vector). A similar law is obtained for the reconstruction of the 
declination as for the amplitude: 



0"dec - 



K dec a y[N 



(21) 



In the case of the Southern site of the PAO, with a dipole in the 
GC direction, one finds an angular power of /fd ec 0.28 rad" 1 , 
or if one prefers K'T} - 200°. Likewise, for the reconstruction 
of the RA, one finds K RA a 0.38 rad" 1 , i.e. K R l A =* 150°. For 
instance, if a = 0.05 (as could be suggested by the AGASA 
data for CRs at ~ 10 18 eV; Hayashida et al. 1999) and for a 
data set of 10 5 events, the angular reconstruction accuracy is 
about 12° for 6 and 9.4° for RA. 



3.3. Comparison with the standard Rayleigh analysis 

In order to compare the power of our method with that of the 
standard analysis in right ascension, we first need to exam- 
ine the link between rih and the dipole amplitude, a. Indeed, 
while the significance of an anisotropy measurement has the 
same form in both cases, Eqs. @ and ( 1201 . it should be kept 
in mind that the first harmonic amplitude in right ascension has 
no clear physical meaning. It would thus be misleading to in- 
terpret Eq. © as demonstrating that the reconstruction power 
of the standard Rayleigh method is 1 / V2 ~ 0.707. To compare 
the reconstruction power of the two methods, we need to evalu- 
ate the significance of the anisotropy measurement/or the same 
value of the dipole amplitude. Therefore, we need to relate rn, 
and a. 

In the standard Rayleigh method, the basic quantities com- 
puted from the data are a and b, given by Eq. Q. These are 
in fact discrete versions of continuous integrals which can be 
identified, according to the relations (|8} and (|9jl, as: 

2 A 2 r 

a — — 2^ cos ifi — — J <l>(u)e(u) cos <pdQ,, (22) 



i=l 



and 



2 A 2 r 

b — — sm Pi - J < I ) (u)e(u) sin cpdCl, (23) 



where N is given by (see Eq.|5}: 

N = J (D(u)e(u)dQ. 



(24) 



The first harmonic amplitude, = Va 2 + b 2 , is thus ob- 
tained by direct integration of the above formulae, over the part 
of the sky where the exposure, S(6) (assumed independent of 



RA), is non zero. For a dipole orientation 84 (in equatorial 
spherical coordinates), one obtains: 



r\h 



C3 a sin #d 



c\ + C2 a cos 0(j 



(25) 



where c\, C2 and C3 are numerical constants which only depend 
on the exposure function of the detector under consideration, 
and are given by: 



ci = 



C2 = 



C3 = 



f 6(0) sin 9 d0 

Je mi „ 

f £(0) sin 61 cos 61 d6 

Je mi „ 

8(0) sin 2 d0 



(26) 



In the case of the PAO detector with the Southern site only, 
we find: a = 0.77, c 2 = -0.32 and c 3 = 0.60. For the full 
Auger detector (both sites), we have: c\ = 1.4, a = 0.028 and 
C3 = 1.1. For a full-sky uniform exposure, the limiting values 
would be c\ — 2, ci — and C3 = tt/2. 

As can be seen, the relation between nti and a (the true, 
physical dipole amplitude) depends on the dipole declination, 
#d = § ~0d, which is precisely not reconstructed by the method. 
It also appears that, as expected, r\h vanishes for<S f / = -n/2 and 
5d = n/2, which indicates that the right ascension analysis is 
totally inadequate to study the anisotropy of the data set in such 
a case. This is obviously because when the dipole vector is ori- 
ented along the rotation axis of the Earth, the flux modulation 
is only in declination. The efficiency of the dipole reconstruc- 
tion with such a method is thus expected to be largest when the 
dipole vector is in the equatorial plane. 

We can now compare the reconstruction powers of the stan- 
dard Rayleigh method and our 3D analysis for different dipole 
orientations. To this purpose, we simply rewrite Eq. @ as a 
function of a: 



n a = -!= rih VA 7 = K lh a Va 7 , 
V2 



with 



^Mh - — — 



1_ rih 
V2 " 



1 



C3 sin6*d 



(27) 



(28) 



C\ + C2 ffCOS f 

Figure |3 shows the dipole amplitude reconstruction power 
both with our method and with the standard one, as a function 
of the dipole declination, in the case of the PAO with one or two 
sites (note that only the Southern site is funded yet, and close 
to completion, while the Northern site in the selection process). 
As can be seen, in the case of a partial sky coverage (Southern 
site only), the first harmonic amplitude method is more power- 
ful for dipole declinations between -60° and +60°. However, it 
should be remembered that only one angle is reconstructed in 
this case, namely the right ascension, and that the inferred is 
not the dipole amplitude (which cannot be derived). For other 
dipole orientations, our method is more powerful (essentially 
because almost no anisotropy remains along RA coordinate), 
and reconstructs entirely the dipole direction, together with its 
true amplitude. 
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Fig. 4. Distributions of the relative error on the reconstructed dipole direction (left), and evolution of the corresponding bias and 
dispersion as a function of the number of events (right), for a dipole of amplitude a - 0.1 oriented towards the GC, observed 
with the Southern site of PAO detector with max = 70°. 
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dipole declination 

Fig. 5. Comparison of the dipole amplitude reconstruction 
powers as a function of the dipole declination, in the case of 
the PAO with one (S) or two sites (N&S). The power of the 
Rayleigh method, Kih, is given by Eq. I28> . 



The reconstruction power, K a , still gets smaller when the 
dipole vector is aligned with the Earth rotation axis, but it does 
not drop to zero, as the 3D analysis also makes use of the decli- 
nation information in the CR arrival directions. The amplitude 
of the declination dependence of K a is limited to a factor of 
~ 2. 

In the case when both PAO sites are available, our recon- 
struction power is always larger than with the Rayleigh method, 
and essentially independent of the dipole orientation. With such 
a full-sky experiment, the reconstructed method presented here 
reduces exactly to that of Sommers (2001), and it is therefore 
natural that we find essentially the same value for K a . The value 



quoted by Sommers (2001) is 0.65, instead of ~ 0.61 in our 
case, but his value was obtained with a smaller data set, and 
as can be seen on Figs. 12) and[3j>, the effective reconstruction 
power for smaller data sets is always slightly higher than the 
asymptotic value obtained at large N (which is the one we plot 
on Fig. 13- Therefore, we can claim excellent agreement with 
previous work when a comparison is possible (i.e. for full sky 
exposure experiments). 

From the experimental point of view, a larger recon- 
struction power implies that a given dipole amplitude will 
be detected earlier, i.e. with a smaller number of events. 
Quantitatively, adding the Northern site to the Pierre Auger 
Observatory will allow one to detect a dipole anisotropy with a 
power between 1 .3 and 2.7 times larger, or a factor of 2 on aver- 
age (for an unknown dipole orientation). According to Eq. ( 1201 . 
this means that its identification will require between 1 .7 and 
7.3 times less events, or a factor of 4 on average. Considering 
that the Northern site will also double the PAO acceptance, and 
thus the rate of events detection, one can expect on average that 
the complete PAO detector will be able to measure a dipole 
anisotropy (at a given significance level) about 8 times quicker 
than just its Southern site (or between 3.4 and 14.6 depending 
on the dipole orientation). 

Concerning the angular accuracy, the reconstruction power 
is plotted in Fig. [(J and[6j> as a function of the dipole declina- 
tion, for either one or two PAO sites. In both cases, we show 
the angular power for the reconstruction of the dipole declina- 
tion, /Tdec an d right ascension, A^a, together with the total an- 
gular power in 3D, K m , The latter is obtained from the 3D an- 
gle, y, between the true and the reconstructed dipole directions, 
given by cos y = cos 6>d cos 6> rec +sin 0d sin 6> rec cos(^d-<A-ec)- F° r 
comparison, we also plotted the pseudo angular power of the 
standard Rayleigh method (relative to the RA projection of the 
dipole direction), which is obtained from the angular accuracy 
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Fig. 6. Angular reconstruction powers of the PAO with one (left) or two sites (right), as a function of the dipole declination. The 
labels refer to the total angular accuracy, and separately to its declination and RA parts, compared to the RA accuracy of the first 
harmonic method. 



on the right ascension recalled in Sect. 12. II cr$ = J2/Nr^ h : 



K RA (lh) = 



< 




1 nil 




V2 a 


l 




C3 sin# c 



y[2 c\ + C2 a cos 6d 

Note that it has the same formal expression as K\h, in 
Eq. d281 i. but it is essentially a different quantity, expressed in 
rad~', while the amplitude reconstruction power is dimension- 
less. 

In the case of the Southern site alone, we first see that the 
reconstruction of the dipole right ascension is always more ac- 
curate with the first harmonic amplitude method, as the full 
statistics is used to derive this parameter alone. On the other 
hand, only the 3D method can reconstruct the dipole declina- 
tion and it is found that the accuracy on this angle is better 
than that on RA for dipole orientations with 6 d ^ -45° or 
5d ~ 45°. As for the dipole amplitude reconstruction, the CR 
flux anisotropy is poorly analysed by the 2D Rayleigh method 
when the dipole vector is at a small angle from the Earth rota- 
tion axis, since only a very small modulation then remains on 
the RA coordinate. 

This transfer of the main flux modulation from one an- 
gular coordinate to the other is also responsible for the anti- 
correlation, observed on Fig. [6J between the reconstruction 
powers for RA and S. The total angular power, given by K tol , 
shows a smaller dependence on the dipole declination, of the 
order of ±15% in the case of the PAO Southern site, with a 
maximum when the dipole is oriented towards one of the poles. 

With the two PAO sites, the angular reconstruction power 
is roughly independent of the dipole declination (as the expo- 
sure is almost flat over the whole sky), and always larger with 
our 3D method than with the standard RA analysis. Even when 
the latter is at its maximum, namely when the dipole vector is 



in the equatorial plane, the 3D reconstruction - i.e. the recon- 
struction of two angles instead of one - is not less accurate than 
the standard 2D Rayleigh reconstruction. Quantitatively, in the 
case of a dipole amplitude of 5% and a data set of 10 5 events, 
the PAO with two sites would give an accuracy of ^ 6.4°, in- 
stead of =s 11° with only one site. 

4. Conclusion 

We have presented a simple numerical method allowing one 
to reconstruct the parameters of a dipolar modulation of an 
otherwise isotropic CR angular distribution. In contrast with 
the standard Rayleigh analysis making use of the distribution 
of CRs in right ascension only, our method enables a full 3D 
reconstruction, providing a direct determination of the under- 
lying dipole right ascension and declination, as well as of its 
true amplitude. We also calculated the relation between the first 
harmonic amplitude in right ascension, r^a, as obtained from 
the 2D Rayleigh method, and the true dipole amplitude, a, as 
a function of the underlying dipole orientation. This relation 
strongly depends on the dipole declination, which cannot be 
reconstructed by the standard method. 

The method presented here is essentially a generalisation 
of that of Sommers (2001), which can now be applied in situ- 
ations where the experimental data do not cover the entire sky. 
In particular, it can be applied with data sets obtained from 
Earth-based experiments detecting CRs arriving from regions 
bounded by lines of equal declination in the equatorial sky. 

In order to quantify the accuracy of the method, we in- 
troduced the concept of reconstruction power, for the dipole 
amplitude as well as for its direction (in 3D space). The sig- 
nificance of an anisotropy measurement can be expressed as 
no- = K a a y/N, where a is the dipole amplitude, N is the num- 
ber of events detected, and K a is the amplitude reconstruction 
power, which depends on the experiment under consideration, 
through its exposure function, £(u), giving the achieved rel- 
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ative exposure of the different regions of the sky. We have 
shown that the reconstruction power depends on the orienta- 
tion of the dipole underlying the CR angular distribution, but 
does not drop to zero when the dipole is aligned with the Earth 
rotation axis, as with the standard 2D Rayleigh method. 

We also investigated the angular accuracy, and showed that 
it can be expressed as cr tot = K tot ayN, where the angular 
reconstruction power, K tot , also depends on the dipole decli- 
nation, but to a smaller extent. In particular, the accuracy of 
the method presented here can be higher than that of the stan- 
dard Rayleigh method, even though it uses the same data set 
to reconstruct one more angular parameter, as well as the true 
(physical) dipole amplitude. Therefore, we believe that this is 
a powerful method to be widely used for angular distribution 
analysis, in particular in CR physics. 

An obvious flaw of the method, however, is that it as- 
sumes a purely dipolar modulation of an isotropic flux (with- 
out any restriction on the dipole amplitude). As was discussed 
in the Introduction, higher order spherical harmonics can be 
present in the underlying angular distribution, which would 
then partly spoil the reconstruction of the dipolar term. In par- 
ticular, a purely quadrupolar distribution could be mistakenly 
reconstructed as a dipole. However, it may first be argued that 
there is a wide range of situations where the dipolar term should 
be dominant in the CR angular distribution, as one might ex- 
pect it to be the first one to be detectable when perfect isotropy 
is lost (unless particular symmetry conditions apply). It should 
be noted also that the situation is the same with the standard 
Rayleigh method, which not only cannot indicate the nature of 
the detected anisotropy (i.e. whether it is a dipole or it involves 
higher order terms), but cannot even provide the dipole ampli- 
tude when the anisotropy is a purely dipolar modulation of an 
isotropic flux. Nevertheless, even if the dominant anisotropy is 
not dipolar, both the first harmonic amplitude analysis in right 
ascension and our 3D method remain useful to determine the 
significance of the observed departure from isotropy. With the 
formalism introduced in this paper, if a "dipole" amplitude a is 
reconstructed from a data set containing N events, the CR an- 
gular distribution can be safely stated to be anisotropic with a 
statistical significance of n a = K a a yfR. Note also that higher 
order spherical harmonics can be investigated through a sys- 
tematic study of the CR angular power spectrum, as reported 
recently by Deligny et al. (2004). 

Finally, as an illustration of our method, we computed 
the angular and amplitude reconstruction powers of the Pierre 
Auger Observatory, comparing the values obtained with the 
Southern site alone or with both Southern and Northern sites. 
We found that the Northern site will lead to a considerable in- 
crease of power, by a factor of ~ 2 on average (depending on 
the underlying dipole orientation) for the reconstruction of both 
the amplitude and the direction (in 3D space). This implies 
that our method will assign a given statistical significance to 
an anisotropy measurement with a data set ~ 4 times smaller if 
two sites are available instead of one. A given measurement 
will thus be reached 8 times quicker with two sites (taking 
into account the corresponding increase in the acceptance) or 
4 times quicker than if the Southern site has its area doubled, 
rather than duplicated in the Northern hemisphere. Conversely, 



for a given size of the overall data set, an anisotropy will be 
measured with typically twice as many sigmas and the angular 
resolution will be twice as more precise with a full sky detector 
than with a partial sky experiment such as the PAO Southern 
site. This may be considered as another important argument in 
favour of the development of the PAO Northern site. 
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